Geospatial Workshop PPSP 2024

Spatial Packages in R

Ensure the following R packages are installed beforehand.

install.packages(c("sf", "terra", "geodata", "rnaturalearth", "spdep",

                   "dplyr", "SpatialEpi", "wbstats", "flexdashboard",

                   "ggplot2", "viridis", "RColorBrewer", "patchwork", "DT",

                   "leaflet", "mapview", "leafpop", "leafsync", "rasterVis"))

install.packages("INLA",
repos = "https://inla.r-inla-download.org/R/stable", dep = TRUE)

Load Packages

For this session, we will load some of the packages needed for plotting basic spatial data

library(sf)

Warning: package 'sf' was built under R version 4.4.1

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

library(mapview)

Warning: package 'mapview' was built under R version 4.4.1

library(janitor)

Warning: package 'janitor' was built under R version 4.4.1

Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

library(readxl)
library(dplyr)

Warning: package 'dplyr' was built under R version 4.4.1

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(ggplot2)

Warning: package 'ggplot2' was built under R version 4.4.1

library(DT)

Warning: package 'DT' was built under R version 4.4.1

Making Maps with R

kelantan <- st_read("kelantan.shp")

Reading layer `kelantan' from data source 
  `C:\Users\ACER\Desktop\geo_workshop\kelantan.shp' using driver `ESRI Shapefile'
Simple feature collection with 66 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 371629.6 ymin: 503028.2 xmax: 519479.6 ymax: 690232.8
Projected CRS: Kertau (RSO) / RSO Malaya (m)

Plot 1

mapview(kelantan[,3])

Plot 2

# Create a custom popup with multiple fields
popup_info <- paste0(
  "<strong>MUKIM: </strong>", kelantan$MUKIM, "<br>",
  "<strong>DAERAH: </strong>", kelantan$DAERAH, "<br>",
  "<strong>Population: </strong>", kelantan$JUM_JANTIN
)
mapView(kelantan, zcol = "JUM_JANTIN", layer.name = "Population", popup = popup_info)

Adding Spatial Attributes

Leptospirosis Data

lepto <- read_xlsx("leptospirosis.xlsx") %>% clean_names()
glimpse(lepto)

Rows: 1,106
Columns: 21
$ diagnosis              <chr> "Leptospirosis", "Leptospirosis", "Leptospirosi…
$ tahun_daftar           <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016,…
$ epid_daftar            <dbl> 6, 9, 11, 12, 18, 18, 19, 17, 21, 23, 34, 34, 3…
$ age                    <dbl> 30, 23, 39, 43, 31, 34, 48, 55, 48, 39, 24, 45,…
$ alamat_semasa_kejadian <chr> "FELDA ARING", "LADANG U&I CIKU 8", "SYARIKAT S…
$ poskod                 <dbl> 18300, 18300, 18300, 18300, 18300, 18300, 18300…
$ latitude_rso           <dbl> 478031, 459494, 441802, 488502, 496760, 441782,…
$ longitude_rso          <dbl> 548141, 564966, 547551, 547701, 539072, 547566,…
$ latitude_wgs           <dbl> 4.944824, 5.034372, 4.897434, 4.945540, 4.93648…
$ longitude_wgs          <dbl> 102.3491, 102.1477, 101.9605, 102.3498, 102.454…
$ notifikasi_status      <chr> "Daftar Kes", "Daftar Kes", "Daftar Kes", "Daft…
$ race                   <chr> "Foreigner", "Foreigner", "Foreigner", "Foreign…
$ kewarganegaraan        <chr> "Bukan Warganegara", "Bukan Warganegara", "Buka…
$ gender                 <chr> "Male", "Male", "Male", "Male", "Male", "Male",…
$ nationality            <chr> "INDONESIA", "INDONESIA", "INDONESIA", "INDONES…
$ klasifikasi_kes        <chr> "Sporadic", "Sporadic", "Sporadic", "Sporadic",…
$ cara_pengesanan_kes    <chr> "Pasif", "Pasif", "Pasif", "Pasif", "Pasif", "P…
$ jenis_rawatan          <chr> "Wad Perubatan", "Jabatan Kecemasan & Kemalanga…
$ daerah                 <chr> "GUA MUSANG", "GUA MUSANG", "GUA MUSANG", "GUA …
$ mukim                  <chr> "CHIKU", "CHIKU", "GALAS", "CHIKU", "CHIKU", "B…
$ negeri                 <chr> "KELANTAN", "KELANTAN", "KELANTAN", "KELANTAN",…

Convert Leptospirosis data to spatial data

Use st_as_sf ( ) function to convert line listing data that contained Lat/Long attributes to spatial object

lepto <- st_as_sf(lepto, 
                    coords = c("longitude_wgs", "latitude_wgs"), 
                    crs = 4326)

Confirm CRS is wgs84

st_crs(lepto)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

Convert shapefile to RSO to match with Kelantan Map

lepto_2 <- st_transform(lepto, 3168)

Plot map to see outlier

ggplot() +
  geom_sf(data = lepto_2) +
  ggtitle("Map of Leptospirosis") +
  theme_bw()

Plot the cases over the Kelantan map

overall_plot <- ggplot() +
  geom_sf(data = kelantan) +   #base map
  geom_sf(data = lepto_2, color = "red", size = 0.5) +  #cases in spatial data
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
  theme_bw() +  
  theme(plot.title = element_text(size = 12),  strip.text = element_text(size = 12)) # cosmetic

overall_plot 

overall_plot + 
  facet_wrap(~tahun_daftar) + #to plot according to year
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5)) +
  theme(plot.title = element_text(size = 12),  strip.text = element_text(size = 12)) +
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022")

Using OSM as basemap

library(leaflet)

Warning: package 'leaflet' was built under R version 4.4.1

leaflet(lepto) |>
  addTiles() |>  # Add default OpenStreetMap tiles
  addCircleMarkers(
    radius = 2,   # Adjust marker size
    color = "red",  # Marker color
    stroke = FALSE, fillOpacity = 0.8)

Aggregated Data

For this exercise, we focus our analysis on Leptospirosis cases reported in 2016 only.

lepto_16 <- lepto_2 %>% filter(tahun_daftar == "2016")

Joint point to polygon

st_crs(lepto_16)

Coordinate Reference System:
  User input: EPSG:3168 
  wkt:
PROJCRS["Kertau (RSO) / RSO Malaya (m)",
    BASEGEOGCRS["Kertau (RSO)",
        DATUM["Kertau (RSO)",
            ELLIPSOID["Everest 1830 (RSO 1969)",6377295.664,300.8017,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4751]],
    CONVERSION["Rectified Skew Orthomorphic Malaya Grid (metre)",
        METHOD["Hotine Oblique Mercator (variant A)",
            ID["EPSG",9812]],
        PARAMETER["Latitude of projection centre",4,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8811]],
        PARAMETER["Longitude of projection centre",102.25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8812]],
        PARAMETER["Azimuth of initial line",323.0257905,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8813]],
        PARAMETER["Angle from Rectified to Skew Grid",323.130102361111,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8814]],
        PARAMETER["Scale factor on initial line",0.99984,
            SCALEUNIT["unity",1],
            ID["EPSG",8815]],
        PARAMETER["False easting",804670.24,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Malaysia - West Malaysia onshore."],
        BBOX[1.21,99.59,6.72,104.6]],
    ID["EPSG",3168]]

Count all leptospirosis cases for each mukim in 2016

count_lepto_sub_yr <- lepto_16 %>% 
  count(daerah, mukim) %>% 
  ungroup()

Joining the count data to base-map of Kelantan

kelantan_count_lepto <- st_join(kelantan, count_lepto_sub_yr)

Map the count of cases according to subdistrict in Kelantan

# Create a custom popup with multiple fields
popup_info <- paste0(
  "<strong>MUKIM: </strong>", kelantan_count_lepto$MUKIM, "<br>",
  "<strong>DAERAH: </strong>", kelantan_count_lepto$DAERAH, "<br>",
  "<strong>Leptospirosis Cases: </strong>", kelantan_count_lepto$n
)

# Display the map with the custom popup
mapView(kelantan_count_lepto, zcol = "n", layer.name = "Count of Leptospirosis cases 2016", popup = popup_info)